Pollinator mediated competition: Comparing traits and environment in reducing phenology overlap of cloud forest plants.
Ben Weinstein
Monday, October 13, 2014
1 Abstract
Cloud forest plant assemblages are among the richness worldwide. Hummingbird visited plants make up a large percentage of the plant assemblage. Previous work has shown that hummingbirds feed on plants that match corollas based on bill length. This suggests that plants with similiar corolla morphologies compete for pollination resources. To minimize competition between species, as well as decrease amount of pollen loss due to heterospecific transfer, there should be selection for minimal overlap in flowering timing. While this has been evaluated within a single genus (Stiles 1975) there has been little work in combining relatedness, environment and morphology in a single framework. My goal is to compare the effect of morphology of competiting species on the abuncance of flowers along a wide elevation gradient.
2 Introduction
2.1 Data
We collected flower abundance for hummingbird visited flowers along a 1300m elevation gradient in northern ecuador.
We measured three flower morphology traits for each species: corolla width, corolla length and height of plant from ground.
2.2 Seasonality
Flower abundance changes throughout the year and along the elevation gradient.
3 Model Formulation
\[ Let \{frac{Abundance}{transect} = \theta\]
Our goal is to explain the the number of flowers per transect based on elevation, date, species, and co-occurring species.
\[ \theta \sim Elevation + Julian Day + Error\]
\[Error \sim \text {Morphology of Co-occurring Species}\]
The goal of the analysis is to measure the effect of this error term on shaping abundance of species flowers. Additional analysis could also look at elevation/distance to a conspecific as promoting abundance.
The abundance of plants is models as a log normal distribution instead of a poisson to have control over the variance term.
\[\theta \sim Poisson(\lambda) \]
3.1 Outline of Steps
This analysis needs to broken down into several tractable steps.
- Model of all plants abundance as a function of intercept + elevation where the model process is considered to follow a poisson.
\[\theta \sim Poisson(\lambda) \]
\[log(\lambda) \sim Intercept + \beta_1 * Elevation \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \]
- Model of all plants abundance as a function of intercept + elevation + Julian Day. The Julian day needs to be transformed?
\[\theta \sim Poisson(\lambda) \]
\[log(\lambda) \sim Intercept + \beta_1 * Elevation + \beta_2 * \text{Julian Day} \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \] \[\beta_2 \sim N(0,.0001) \]
- Model of all plants abundance as a function of intercept + elevation + Julian Day. Sigma is modeled as function of Morphology of co-occurring species. This now considered a hierichical model.
\[\theta \sim LogN(\mu,\sigma) \]
\[\mu \sim Intercept + \beta_1 * Elevation + \beta_2 * \text{Julian Day}. \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \] \[ \text{Julian Day} \sim U(0,365) \] \[\beta_2 \sim N(0,.0001) \]
With the errors: \[\sigma \sim Gamma(lambda,.001)\] \[ \lambda \sim Intercept + \beta_3 * \text{Morphology of Co-occurring Species}\] \[ \beta_3 \sim N(0,0.0001)\]
- Model of all plants abundance as a function of intercept + elevation + Julian Day. Sigma is modeled as function of Morphology of co-occurring species. A random effect of species is added for each species. Is this written:
\[\theta \sim LogN(\mu,\sigma) \]
\[\mu \sim Intercept + \beta_1 * Elevation + \beta_2 * \text{Julian Day}. \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \] \[ \text{Julian Day} \sim U(0,365) \] \[\beta_2 \sim N(0,.0001) \]
With the errors: \[\sigma \sim Gamma(lambda,.001)\] \[ \lambda \sim Intercept + \beta_3 * \text{Morphology of Co-occurring Species_ij + Species_i}\] \[ \beta_3 \sim N(0,0.0001)\]
Where the random effect of Species is normally distribution \[Species_i \sim N(0,0.0001)\]
3.2 Analysis
d<-read.csv("C:/Users/Ben/Dropbox/Thesis/Maquipucuna_SantaLucia/Results/FlowerTransects/FlowerTransectClean.csv",row.names=1)
#remove missing or infinite data
d<-d[is.finite(d$Total_Flowers),]
d<-d[is.finite(d$ele),]
d<-d[!d$Total_Flowers==0,]
#Get the transect totals
s.dat<-split(dat,list(dat$ID),drop=TRUE)
total_flower<-rbind.fill(lapply(s.dat,function(x){
#mean elevation
Elev<-mean(as.numeric(strsplit(as.character(x$Transect_R),"_")[[1]]))
#Date
Date<-unique(x$Date)
#Total Number of Flowers
fl<-round(sum(x$Total_Flowers,na.rm=TRUE))
data.frame(Elev,Date,Flowers=fl)
}))
#remove NA
total_flower<-total_flower[!is.na(total_flower$Flowers),]3.3 Step 1 Bayesian model of elevation of all species combined using a Polynomial Poisson
3.3.1 Visualize data
ggplot(total_flower,aes(x=Elev,y=Flowers)) + geom_point() + geom_smooth(method="glm",family="poisson",formula=y~poly(x,2)) + ggtitle("Polynomial Poisson GLM")
#############################
# Load all the libraries
#############################
library(R2jags)
setwd("C:/Users/Ben/Documents/PhenoBayes/")
source("Model1.R")
# cat("
# model {
# for(i in 1:length(y)){
# y[i] ~ dpois(mu[i])
# mu[i] <- alpha + beta * Elevation[i] + beta2*pow(Elevation[i],2)
# }
# alpha ~ dnorm(0,0.001)
# beta ~ dnorm(0,0.001)
# beta2 ~ dnorm(0,0.001)
#
# }",fill = TRUE)
#
# sink()
#############################
#Make a list where you include all the data the model will need to run
#############################
Dat <- list(
y=total_flower$Flowers,
Elevation=total_flower$Elev
)
#############################
#Initial Values
#############################
#Make a column vector with the names of the parameters you want to track
#############################
ParsStage <- c("alpha","beta",'beta2')
#############################
#Set the variables for the MCMC
#############################
ni <- 5000 # number of draws from the posterior
nt <- 1 #thinning rate
nb <- 0 # number to discard for burn-in
nc <- 2 # number of chains
#############################
#Run the jags function to run the code
#############################
m = jags(
n.chains=nc,
model.file="flower.jags",
working.directory=getwd(),
data=Dat,
parameters.to.save=ParsStage,
n.thin=nt,
n.iter=ni,
n.burnin=nb,
DIC=T)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 344
##
## Initializing model3.3.2 View model
m## Inference for Bugs model at "flower.jags", fit using jags,
## 2 chains, each with 5000 iterations (first 0 discarded)
## n.sims = 10000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75%
## alpha -20.616 0.126 -20.858 -2.07e+01 -20.615 -20.533
## beta 0.031 0.000 0.031 3.10e-02 0.031 0.031
## beta2 0.000 0.000 0.000 0.00e+00 0.000 0.000
## deviance 254443.735 724.961 253036.804 2.54e+05 254455.000 254918.071
## 97.5% Rhat n.eff
## alpha -20.371 1.001 10000
## beta 0.031 1.001 10000
## beta2 0.000 1.000 1
## deviance 255842.182 1.001 10000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 262798.4 and DIC = 517242.2
## DIC is an estimate of expected predictive error (lower deviance is better).3.3.3 Reshape data
require(reshape2)
pars<-melt(m$BUGSoutput$sims.array[-c(1:2000),,c('alpha','beta','beta2')])
colnames(pars)<-c("Draw","Chain","parameter","estimate")3.4 Visualize chains
require(ggplot2)
ggplot(pars,aes(x=Draw,y=estimate,col=factor(Chain))) + geom_point(size=.5) + geom_line() + ggtitle("Chains") + facet_wrap(~parameter,scales="free") + theme_bw() + labs(col="Chain")
3.5 Histogram of Posteriors
require(ggplot2)
ggplot(pars,aes(x=estimate)) + geom_histogram() + ggtitle("Estimate of parameters") + facet_wrap(~parameter,scale="free") + theme_bw() 
3.6 Plot the data and the bayesian best fit line
#reformat for to sample posteriors
post<-dcast(pars,...~parameter)
#sample 1000 rows from the posterior
sample.post<-post[sample(1:nrow(post),1000),]
head(sample.post)## Draw Chain alpha beta beta2
## 1049 525 1 -20.42 0.03082 -8.313e-06
## 4501 2251 1 -20.42 0.03083 -8.317e-06
## 3877 1939 1 -20.62 0.03103 -8.365e-06
## 3480 1740 2 -20.56 0.03095 -8.342e-06
## 276 138 2 -20.66 0.03107 -8.378e-06
## 2293 1147 1 -20.55 0.03095 -8.345e-06#for each year multiply the value by the sample of the posterior
confidence.years<-lapply(total_flower$Elev,function(x){
exp(sample.post$alpha + sample.post$beta * x + sample.post$beta2*x^2)
})
names(confidence.years)<-total_flower$Elev
means<-sapply(confidence.years,mean)
#calculate the central 95th CI interval
quantile_elev<-t(sapply(confidence.years,function(x){
quantile(x,c(0.025,0.975))}
))
#reformat for ggplot
confidence<-melt(confidence.years)
#name and make year a number
colnames(confidence)<-c("Estimate","Elev")
confidence$Elev<-as.numeric(confidence$Elev)
my<-data.frame(total_flower$Elev,means,quantile_elev)
colnames(my)<-c("Elev","mean","lower","upper")
p<-ggplot(confidence,aes(x=Elev,y=Estimate))+ theme_bw() + geom_line(data=my,aes(x=Elev,y=mean,group=1),col="red") + geom_line(data=my,aes(x=Elev,y=lower,group=1),col="blue",linetype='dashed') + geom_line(data=my,aes(x=Elev,y=upper,group=1),col="blue",linetype='dashed') + ggtitle("Confidence Interval for y~Year")
#Add original data
p <-p + geom_point(data=total_flower,aes(x=Elev,y=Flowers),size=3)
p